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The threshold probability for successful error correction of a quantum system is not a single, 
absolute value. Not only is it contingent on the noise model, but what the noise model is believed 
to be. This is exhibited for the 2D Toric code subject to local noise where, via a mapping to the 
2D random bond Ising model, suitable thresholds can be established. The analytic calculation is 
founded on an assumption which is known to give a good approximate result. We also find the first 
order correction to this result based on renormalisation group calculations, and numerically evaluate 
lower bounds based on explicitly implementing an error correction strategy. 



The possibility of error correcting quantum systems 
constitutes a cornerstone of quantum computation; with- 
out it, information processing in the quantum regime 
would be inherently limited to small systems, providing 
minimal computational advantage. Surface codes, specif- 
ically the Toric code in 2D pQ , provided another massive 
step forward, since error correction can be performed 
fault tolerantly with a finite threshold [2J, i.e. there is 
a non-zero per-qubit error rate below which error correc- 
tion is successful even in the presence of faulty measure- 
ments, with a spatially local structure. This removes the 
need for the complicated hierarchy of concatenated codes 
present in the original threshold derivations which, while 
technically impressive, were highly impractical. 

The simplicity, topology and locality of the Toric code 
have also engendered the idea of a quantum memory [3]-[5] 
wherein a quantum state is encoded in the ground state 
of a local Hamiltonian, and the structure of that Hamilto- 
nian prevents any locally occurring errors from aggregat- 
ing into topologically non-trivial errors which could cor- 
rupt the encoding, thereby minimising the rate at which 
error correction must be applied. Errors in the Hamil- 
tonian, or the evolution of external errors driven by the 
Hamiltonian, typically do not manifest as independent lo- 
cal rotations. Nevertheless, they must be corrected even 
if their distribution is not explicitly known. 

In this paper, we derive a series of new thresholds for 
the Toric code in 2D which elucidate the influence of the 
difference between the true noise model and our assump- 
tions about its nature. For simplicity of exposition, error 
correction is assumed to be perfect, although imperfect 
measurements can be incorporated in principle. 

The main technical tool is the previously established 
connection between syndrome measurements on a noisy 
Toric code and the random bond Ising model (RBIM) in 
2D [2J. The phase transition of the RBIM locates the 
critical threshold of the Toric code. In [2] , this connec- 
tion was established for a noise model parametrised only 
by p, the error rate, and assumed that the noise model 
was known. We extend this to a noise model that con- 
tains two error rates px and pz, and our assumptions 
about what these values are, px and pz- The critical 



region of the RBIM can be determined by an ansatz [5] 
and improved upon by a renormalisation style expansion 
[TJ. Numerical verification comes via the simulation of 
the performance of an explicit error correction procedure 
(thereby yielding a lower bound on the threshold), mini- 
mum weight perfect matching. 

Toric Code: The Toric code pQ is the quintessential 
example of a planar code; one which can be defined for 
arbitrarily many qubits with a correspondingly increas- 
ing code distance, and which has a finite fault-tolerant 
threshold without requiring the additional complications 
of concatenation. We consider here, as in [3], a rotated 
version akin to the Wen code |S]. This only changes the 
interpretation of the types of noise that affect each site, 
making them more natural. Starting from an N x N 
square lattice with periodic boundary conditions, place a 
qubit in the middle of each edge. For each vertex v and 
face / of the lattice, there are four neighbouring qubits, 
two on horizontal edges Eh and two on vertical edges 
Ey . The corresponding stabilizers are defined as 

K v = J] Z e Y[ X e , K f = ]J X e l[ Z e . 

All terms mutually commute and each has eigenvalues 
±1. The space of Toric code states l^y) for i,j £ 
{0,1} are defined by the relations K v = and 
Kf \k>ij) = for all v, f. There are 2N 2 qubits and 

2A 2 — 2 independent stabilizers, leaving a four-fold degen- 
eracy (the indices that represents two logical qubits, 
for which the logical Pauli operators correspond to long- 
range operations along a single row or column, looping 
around the entire torus. The logical Z operators are two 
inequivalent loops of Paulis along a row; one made up 
of X rotations along a set of horizontal edges and the 
other Z rotations along a set of vertical edges. The logi- 
cal X operators similarly utilise rotations on two different 
columns. Starting from a logical state l^y), and applying 
continuous segments of X and Z operators, it is possible 
to form closed loops (meaning that all stabilizers return 
+1 expectation). Provided those loops are topologically 
trivial (i.e. they don't form loops around the torus), the 
state is the same as the original one while non-trivial 
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loops correspond to logical errors. This degeneracy of 
the code means that if a large set of errors has arisen, 
it is not necessary to establish from the locations of the 
anyons exactly which errors occurred in order to correct 
for them; one only has to form the closed loops which are 
most likely to be trivial. 

Henceforth, we shall restrict the errors to being such 
that X and Z errors act independently on each site, with 
probabilities px and pz respectively. This assumption, 
which is relaxed in the Appendix, means that there are 
two independent sets of anyons that can never interact; 
those detected by the {Kf} and {K v } respectively. It is 
sufficient to consider just one of these sets, as the other 
will be the same by symmetry. We concentrate on the 
{K v } centred on the vertices. 

Statistical Mechanics: In [2J, a connection was proven 
between the ability to correct errors arising on the Toric 
code and the existence of a phase transition in the ran- 
dom bond Ising model. This was done for both perfect 
(the 2D RBIM) and imperfect stabilizer measurement 
(3D RBIM), assuming that X and Z errors each occur 
with probability p independently on each lattice site. We 
will now extend this argument for the case of perfect sta- 
bilizer measurement. To emphasise the difference with 
the actual error rates {px and pz), the assumed error 
rates are denoted by px and pz respectively. This is an 
important distinction because one should only be able to 
achieve the optimal recovery specified by [SJ if the na- 
ture of the noise is known exactly, but this cannot be 
extracted from the stabilizer measurements; only some 
net error rate p is apparent: 

(l-2p) 2 = (l-2p x )(l-2pz). 

Having extracted the error syndrome by measuring the 
stabilizers, there are only four inequivalent consequences 
(1, Xf, Z% and X^Z%) of the correction. Four possible 
corrections, error strings Ei, i = . . . 3, can be identified 
which would give topologically inequivalent results and 
need to be assigned a likelihood of having arisen accord- 
ing to the assumptions on the noise model 

Pi = ^p(E,UC), 

ces 

where the set S corresponds to all trivial loops. If error 
correction is possible, then the expectation of the prob- 
ability of getting the right answer over all actual error 
configurations should tend to 1 in the limit of large sys- 
tem size, N, while the other probabilities should vanish. 
When error correction fails, all the pi will be similar. 

Let r° be a set of variables ±1 for each qubit corre- 
sponding to whether or not a rotation is applied in the 
correction E n . Similarly, Tq is the set due to fijUC. The 
set of all closed loops is conveniently described by intro- 
ducing variables cr, € {±1} for each vertex of the dual 
lattice [21 IHj ■ A qubit q on an edge of the primal lattice 
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FIG. 1: Phase diagram of the 2D RBIM. The dashed line 
separates the ferromagnetic (FM) and paramagnetic (PM) 
phases. Where the Nishimori line intersects this is the multi- 
critical point. 

has two neighbouring vertices of the dual lattice, v\ and 
v 2- T c q = T q a v\ a v\ i where it now suffices to sum over 
the variables <7j without restriction. For probabilities 

_ J Pz qeV 
Pq {PX qeH 

the probability of a given error string is 

KrS)=II( 1 -^) (1+TS - 9)/2 ^ 1 " T " 5)/2 - 
Removing a common factor and defining 

1 -PX = e 2J H ^-PZ = e 2J v 

Px Pz 
then the probability po is proportional to 

Zo = J2e H ^ 

with 

H{°) = J2( T q J H)v qi e q2 + J2(T°J v )cr qi a q2 . 

This is the Hamiltonian of the ±J random bond Ising 
model on a square lattice, where the vector r° arises 
from the actual errors that occurred, and the coupling 
strengths constitute our assignment of the likelihood of 
different configurations. The transition in behaviour of 
the probabilities between successful correction (in asymp- 
totically all cases of instances of the syndromes) and fail- 
ure corresponds to a discontinuity of the free energy of 
this model F = In Z where Z = J2i This is the well 
studied phase transition in the 2D RBIM. The salient 
details of this study are reproduced in the Appendix. 

The most studied instance, px — Pz = Px = Pz, 
transforms to the Nishimori line of the RBIM, and reveals 
the critical probability quoted in [2] via correspondence 
to the multicritical point of the RBIM. 

RBIM: The analysis conveys an interest in the crit- 
ical point of the RBIM, not only along the Nishimori 



3 



line, and for asymmetries between horizontal and verti- 
cal bond strengths. The discontinuity in the free energy 
is located by appealing to the replica method: 



F = In Z = lim 

n->0 



(Z n ) - 1 



The angular brackets denote the average over possible 
configurations of the r° coefficients. 

In the absence of randomness, the 2D Ising model is 
self dual, Z(J Hl J v ) = Z*(J H ,J V ) = Z(J H ,J V ) (the 
dual is denoted by *), where 



-2JT 



tanh( Jy) 



-2J* 



tanh( Jh )■ 



These relations convey that under the duality transfor- 
mation [S], horizontal edges are mapped to vertical ones 
and vice versa. In parameter regimes where the model 
exhibits only one phase transition, the duality relation 
reveals that point, Jh = Jh and Jv — Jv (here, one 
relation implies the other). It is appealing to expect 
that the random model should exhibit similar properties 
and, indeed, the n replicas are self dual but with a larger 
parameter set [6J: {x k ,x k } k=0 for the primal partition 
function and {x* k H , x* k v }" =0 for the dual. 



p x e njH + (l-p x )e 
2 n/2 costi 1 (J H ) 



-nJ H 



(1) 

(2) 



The V versions are equivalent, with the replacements 
Px ^ Pz and Jh h> Jv ■ Unfortunately, there is no sat 
isfying assignment such that x k = x* k v and x\ = X n 
for all k. Instead, it is traditional to proceed with the 
conjecture that the critical point is given by 



H 



Taking the limit n 



h v *h *v 
yields 



(3) 



PH log 2 PH + (1 -Pi?)l0g 2 (l ~Ph) 

+Pv log 2 Py + (1 -pv)log 2 (l ~Pv) = -1 (4) 

The above conjecture has been specifically postulated 
to work at the multicritical points of the RBIM (px = px 
and pz — Pz)-, being well justified by the existence of var- 
ious symmetries, and well tested numerically, also in the 
asymmetric case [BJ. Outside this regime, the results are 
less trusted. In the homogeneous case (pz = Px = P and 
Pz = Px = p), detailed studies [7] justified that above 
the Nishimori line (p > p, see Fig. [T]), the approximation 
was a good one. However, below the Nishimori line, the 
results were inaccurate. On the other hand, the critical 
probability on the Nishimori line gives an upper bound 
to this critical probability (and is quite tight as the phase 
boundary is roughly vertical). In the asymmetric case, 
there is not even a clear notion of being above the Nishi- 
mori line as the working regime 

(1 - 2p x ) 2 (l - 2p z ) 2 = (1 - 2 Px ) 2 (l - 2 PZ ) 2 
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FIG. 2: Comparison for error threshold when the error model 
is known (px = Px,pz = Pz) and when errors are assumed 
equally likely, px = Pz- Plotted for comparison are numeri- 
cal thresholds from minimum weight perfect matching simu- 
lations with a lattice size iV = 100. 

has px > Px and pz < pz or vice versa. However, this 
quandary can be resolved in two relevant cases. 

Firstly, the absolute best that any error correction 
strategy could possibly perform is given precisely by the 
multicritical point, corresponding to perfect knowledge 
of the error model; px = px and pz = pz- 

With the exception of px = [TS], the measurement 
statistics only reveal the frequency p of stabilizers be- 
ing — 1. Without any further knowledge, one should as- 
sume that px = pz- How detrimental is this to the 
threshold? Under this condition, Eq. Q reduces to an 
effective homogeneous system with 2p = px + pz and 
(1 - 2p) 2 = (1 - 2p x )(l - 2p z ). Aside from p x = pz, 
the point lies above the Nishimori line, and the conjec- 
ture of Eq. Q is still reliable [7] . Alternatively, having 
reduced to the homogeneous case, the critical probabil- 
ity cannot be larger than that at the multicritical point, 
Pc ~ 0.1092. Hence, setting p — p yields px +pz < 2pc- 

Although the conjecture of Eq. (J3| compares 
favourably with numerical estimates on the square lat- 
tice, there is known to be a discrepancy with some exact 
renormalisation group calculations on hierarchical lat- 
tices (10) . To account for this, a renormalisation inspired 
expansion has been introduced to account for corrections. 
While most natural for hierarchical lattices [11], it has 
been extended to square lattices and achieves an even 
tighter match with previous numerical results using only 
a first order correction [?]• This technique proceeds by 
replacing the term e J " in Xq (and, similarly, the term 



-Jh 



in Xq ) with an equivalent effective weight 
arising from an averaging effect over several neighbour- 
ing spins. This averaging involves horizontal and vertical 
edges equally, thereby simplifying the conjectured posi- 
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tion of the critical point to xq = Xq. The larger the block 
that is averaged over, the better the anticipated approx- 
imation of the position of the critical point. The first 
order approximation specified by Fig. 4 of [7 locates the 
critical point at 

2 = \ £ (1 + ^(1- 2p) 4 )log 2 (l + ?7(1 -2p) 4 ) 
n=±i 

1 2 /2\ 

_ 2jX/( ) a nr, l {PX,Pz)t0g 2 a nm (px,Pz) 



n— m— 



where 



(1 - r)"- 2 (l - s) r 



(1 - r) n (l - s) r 

r n— 2 c m—2 



The threshold values are barely impacted, although Fig. 
[2] indicates that they can surpass the Gilbert- Varshamov 
bound for stabilizer codes (which coincides with Eq. Q). 
As the code is degenerate, there is no contradiction [12] . 

Minimum Weight Perfect Matching: Given the near- 
vertical phase boundary of the RBIM below the Nishi- 
mori line, any error correction strategy which assumes a 
lower 'temperature' has an almost identical critical prob- 
ability. In particular, the temperature case corresponds 
to correcting by minimum weight perfect matching [20] : 
an efficient algorithm which is readily implemented (how- 
ever, it cannot deal with correlations between the proba- 
bilities of X and Z errors, and a polynomial scaling may 
be insufficient [E]). As such, this algorithm gives a lower 
bound on the threshold fidelities for comparison, and is 
anticipated to be quite tight. The distance measure is 
the log of the probability of a given error path, i.e. a sum 
over pair-wise matching of anyons separated by In and 
ly. Hence, each pair of vertices is assigned a weight 



1 H 



In 



Px 



1 ~Px 



lv In 



Pz 



1-Pz 



and the total weight is minimised [2J. The points in Fig. 
[2] represent a Monte Carlo simulation in a system with 
N = 100 based on the Blossom V algorithm (T3]. For a 
fixed ratio px/Pz, the fraction of 500 runs giving a logi- 
cal error was computed for varying error rates, enabling 
determination of the failure probability (the threshold at 
which a transition in logical error rate from to 50% oc- 
curs). Similar numerics, for a perfectly identified error 
model, are present in [15] . 

Conclusions: In this paper, we have demonstrated how 
any discrepancy between the actual and assumed noise 
models affects the error correcting threshold of the Toric 
code in 2D. For perfectly identified errors, the rotated 
Toric code that we examined (approximately) saturates, 
and even surpasses, the hashing bound, for all asymme- 
tries, without resorting to random lattices [TS]. In com- 
parison, the non-rotated version has the much weaker 
threshold of max(px,pz) = Pc- The present results are 



particularly important in the memory setting. In [3] it 
was demonstrated that a randomised Toric code could 
have a polynomial (in N) survival time between rounds 
of error correction. Any lack of knowledge obliterates 
that polynomial survival time, making it only a constant, 
due to the shifted threshold value. As detailed in the Ap- 
pendix, correlations due to Y errors can be incorporated 
into the analysis [16] , with the equivalent of Eq. Q be- 
ing given in Eq. (A.l). In principle, the present results 



can be generalised to the situation of noisy measurements 
[HE], adding a third dimension to model. 

Acknowledgements: This work was supported by the 
National Research Foundation & Ministry of Education, 
Singapore. 



[1] A. Kitaev, Annals of Physics 321, 2 (2006). 

[2] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, J. 
Math. Phys. 43, 4452 (2002). 

[3] A. Kay, Phys. Rev. Lett. 107, 270502 (2011). 

[4] F. Pastawski, A. Kay, N. Schuch, and I. Cirac, Phys. 
Rev. Lett. 103, 080501 (2009); F. Pastawski, A. Kay, 
N. Schuch, and I. Cirac, Quantum Inf. Comput. 10, 580 

(2010) ; A. Kay, Phys. Rev. Lett. 102, 070503 (2009). 

[5] S. Bravyi and J. Haah, Phys. Rev. Lett. 107, 150504 

(2011) ; S. Bravyi, M. B. Hastings, and S. Michalakis, J. 
Math. Phys. 51, 093512 (2010); J. R. Wootton and J. K. 
Pachos, Phys. Rev. Lett. 107, 030503 (2011); C. Stark, 
L. Pollet, A. Imamoglu, and R. Renner, Phys. Rev. Lett. 
107, 030504 (2011). 

[6] K. Takeda, T. Sasamoto, and H. Nishimori, J. Phys. A: 

Math. Gen. 38, 3751 (2005). 
[7] M. Ohzeki, Phys. Rev. E 79, 021129 (2009). 
[8] X.-G. Wen, Phys. Rev. Lett. 90, 016803 (2003). 
[9] F. Y. Wu and Y. K. Wang, J. Math. Phys. 17, 439 (1976). 
[10] M. Hinczewski and A. N. Berker, Phys. Rev. B 72, 
144402 (2005). 

[11] M. Ohzeki, H. Nishimori, and A. N. Berker, Phys. Rev. 

E 77, 061116 (2008). 
[12] D. P. DiVincenzo, P. W. Shor, and J. A. Smolin, Phys. 

Rev. A 57, 830839 (1998). 
[13] G. Duclos-Cianci and D. Poulin, Information Theory 

Workshop, IEEE , 1-5 (2010). 
[14] V. Kolmogorov, Math. Prog. Comp. 1, 4367 (2009). 
[15] B. Rothlisberger, J. R. Wootton, R. M. Heath, J. K. 

Pachos, and D. Loss, Phys. Rev. A 85, 022313 (2012). 
[16] H. Bombin, R. S. Andrist, M. Ohzeki, H. G. Katzgraber, 

and M. A. Martin-Delgado, Phys. Rev. X 2, 021004 

(2012) . 

[17] K. Takeda and H. Nishimori, Nuclear Physics B 686, 
377396 (2004). 

[18] G. Duclos-Cianci and D. Poulin, Phys. Rev. Lett. 104, 
050504 (2010). 

[19] Every column has an even number of —1 measurement 
results of the stabilizers. The chance of that happening 
if px 7^ is vanishingly small. 

[20] Zero temperature corresponds to the limit J — > 00, mean- 
ing the only term in Zo that is worth considering is the 
<t that gives the smallest value. 



5 



Appendix: Generalised Noise Model 

As a generalisation to the main body of the paper, we want to derive the error correcting threshold of the Toric 
code in 2D when subject to a per-qubit noise model of 

£{p) = (l-qx-qy- qz)p + qxXpX + qyYpY + q z ZpZ 

when we assume that the corresponding probabilities are qx, qy and qz- For comparison, the previous distribution 
had 

qx=px(i-pz) qz=pz(l-px) qv=PxPz- 

With Y errors present, it will not be possible to divide the system into two independent systems. Let us denote by L 
the lattice on which the toric code is defined (with the qubits in the middles of the edges) . Let L\ be a copy of L and 
1/2 be the dual of L\ (i.e. in the case of the periodic square lattice, the same lattice but shifted both horizontally and 
vertically by half a unit). With the edges q (corresponding to a qubit on L) of each of the lattices Li we associate a 
variable T % q £ ±1 in the following way: 

• If q is a vertical bond of L, an X error on q is represented by = — 1. 

• If q is a vertical bond of L, a Z error on q is represented by = — 1. 

• If q is a horizontal bond of L, an X error on q is represented by = — 1. 

• If q is a horizontal bond of L, a Z error on q is represented by = — 1. 

Note that this means that X errors are specified by horizontal edges of L\ and L2, while Z errors are specified by 
vertical edges. For a given error configuration r = (r , r ), we're interested in assessing our belief of it having arisen 
from a given sector, i.e. 

Po = p{tUt Ci Ut C2 ) 

T C ieS 1 ,T C 2£S 2 

where Si and £2 are closed loops around L\ and L2 respectively. This p is proportional to 
where the ' denotes the sum being restricted to closed loops and 

e Aj H = (i-qx - qy - qz)qz ^ Jv = (l-gx-qy- qz)gx ^ Jy = {l-gx - qy - qz)qy 
qyqx qyqz qxqz 

A convenient way to parametrise the sum over closed loops is to instead let an edge q specify the vertices i and j of 
the corresponding edge of L\ and i' and f of Li (see Fig.^J), and to write t^ 1 = (TiCTj, t^ 2 — Oi'Oy . When one takes 
the product of values 1 around a closed loop, every vertex appears an even number of times, and since af = 1, this 
gives the correct product. Hence, we have 

Z Q = Ye H 

a 

with 

H = T qJH&iVj + TgJvCi'G j> +TqTqJy(Tiaj(Ti>aji. 

The T l q have a sign distribution specified by the error model. This is exactly the model derived in |16j . and it 
corresponds to two planes of Ising models coupled by some 4-body terms. It is equivalent to the 8-vertex model. 
We are now interested in determining its phase transitions as this will convey the limit of when error correction is 
possible. In order to do this, the first thing to establish is duality of the non-random version. This can then be 
extended to the replica method in order to account for the configurational average. We note that if we take the 
previous parametrisation of the probabilities, Jy = 0, we are returned to having two uncoupled 2D Ising lattices. 



FIG. 3: Each qubit, q lies on an edge of both the primal (solid) and dual (dashed) lattices. The corresponding vertex pairs 
that comprise the edge are and respectively. The extended model consists of ±1 variables on each vertex of the 

two lattices, with pairwise interactions between and pairs and a four body interaction between 

Duality 

Duality is an important concept for us because if a model's partition function Z(J) for a particular set of parameters 
can be mapped on to itself for a different set of parameters Z(J*), then along a path in parameter space that contains 
its duals, if we know that there is a single discontinuity in Z, it can only happen at J — J* (otherwise discontinuities 
would have to exist in symmetric pairs about this point). For simple noise models in which a single parameter 
describes the strength of the noise, as the strength increases the system must pass through a critical point demarking 
the failure of error correction, and there will only be one such point - adding extra noise should never bring us round 
full circle from correction being impossible to suddenly being able to correct again! 

Let us approach the task of proving duality in the non-random model by expressing 

z = E' IK 

t c i,t c 2 q 

This is just writing exactly what we had before but with 

U q = e * " 1 +T q JVT t + T q T q J Y T q T q , 

Defining the dual functions as 

2u* q = e T l JH+r « Jv+ < T l JY + (-l) s2 q e T q J H-^Jv-ry q J Y + (_ 1 )«J e -T, 1 Ji r +r,Vv-r, 1 Tf Jy + ^yl+s'^-r^ Jn-r 2 q Jv+rtf J Y ^ 

for sj, s q £ {0, 1}, then u q can be expressed as 

»J,«?e{±i} 

having replaced Oi with 2e^ — 1 such that e\ E {0, 1}. Consider what happens when we examine the product of u q , 
expressed as above, for all edges of L leading out of a given vertex i, and performing the sum in Z over that particular 
£i. The terms are unless 

s l mod 2 = 0, 

3 

which is entirely equivalent to the product of s q having to be +1 around vertices or, in other words, around closed 
loops of the dual lattice. Hence, 

*=E'IK 

s 1 ^ 2 q 

meaning that the model is self dual. For instance, the values 

^Jh+Jv+Jy & Jh~Jv-Jy 6 —Jh+Jv-Jy e —JH-Jv+Jy\ 

(11 i 1 \ 

1 1-11-1 

2 11-1-1 " 
\1 -1 -1 1 J 
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are equal if Jh = Jv and e 2Jy = sinh(2Jff). This is the only instance in which the model is self-dual. However, it 
also exhibits a pseudo mutual duality between parameter sets 

(Jh, Jv, Jy) (Jv, Jh, Jy) 

and all other permutations. This is the equivalent of the mutual duality applied to the anisotropic RBIM in the body 
of the paper and [B], and is succinctly expressed as the contour 

qx + qy + qz = \- 

This symmetry property helps to justify the applicability of the assumption we will make to determine the critical 
point of the random model [6J . 

Replica Method 

The above duality was proven without any randomness present. In order to deal with the randomness of the bonds, 
we must take a configurational average, and use it to determine any discontinuity in the free energy, In Z . For n 
copies, Z n = Z® n is dual to with a corresponding set of parameters w®" and u f ® n . If, after averaging over the 
possible values of r*, these two vectors are equal, there is a self-dual point. Unfortunately, this does not happen 
simultaneously for all elements of the vector. Instead, we make the zero order assumption that the critical point 
occurs when the first element of each vector (after averaging) is equal. 

lim x n = x* n . 

n— >0 

Note that changing the values of the permutes the elements of the vector u® n , which means that the possible first 
elements are the n th power of the elements of u. Hence, 

x n = 0--Qx~9r- q~z,qz, qx,qy) ■ u n - 

On the other hand, the first element of the vector u*® n does not change as values are randomly selected. Hence, 

x * _ J_ ( e Jir+Jv+JY i £ Jh-Jv-Jy i 6 ~Jh+Jv~Jy i e -JH-Jv+Jy\ n 
n on V / 

Expanding this, and taking the limit, gives the equivalent of Eq. Q for the generalised model: 

{l-qx-Qr- qz,qz,qx,qv) ■ log 2 (l -qx-Qr- qz,qz,qx,qr) = -l, (A.i) 

which coincides with the Hashing bound for stabilizer codes subject to Pauli error channels. This can be anticipated 
to be a good approximation whenever the error rates have been perfectly identified. 

Now we can ask a variety of questions. For instance, if qx = qy = qz = qx = qy — qz, i.e. depolarising noise that 
we have perfectly identified, one has to solve the equation 

(1 - 3g)log 2 (l - 3q) + 3q\og 2 {q) = -1 

to find iq = 0.18929, which replicates the value given in [16] , where this value was also numerically verified. 

How effective is minimum weight perfect matching at correcting depolarising errors? Since it is not capable of taking 
the correlations introduced by Y into account, it is not expected to be tight with the optimal correction. However, a 
good estimate is to set qy = p 2 , qx = qz = p(l — p), qx — <jy = qz = q, and find the critical point. It's described by 

-\ = (l-2«)log 2 (l-p) + 2 ? log 2 (p), 

and the largest value of q is given by setting p = 2q, which reveals that 3g = 0.165. This compares favourably with 
previous results [TB] . 

One could improve the accuracy of these conclusions by utilising the higher order expansions [7]. However, at low 
orders (i.e. those which are computationally accessible), the differences are negligible, as evidenced in Fig. [2] 



